Getting Started

Assuming we have some extra time, this additional exercise is going to explore the species location data in more detail. We will produce two interactive map examples that you can open in your web browser. These types of simple web maps and applications allow you to share your data quickly and easily with project partners, colleagues, and the general public, with a lot of control over style and content. These examples are also meant to emphasize the amazing capabilities of existing APIs that you can leverage to access high resolution imagery and data visualization functionality that would otherwise be cost prohibitive to purchase or code from scratch.

Load packages

We’ve already loaded a few of these packages in the previous example, but let’s include everything we need for this example just to be safe.

# Packages for Leaflet Mapping Example
library(rgdal)
library(raster)
library(sp)
library(sf)
library(dplyr)
library(tidyr)
library(ggplot2)
library(leaflet)
library(leaflet.extras)
library(leaflet.opacity)
library(htmltools)
library(htmlwidgets)

# Packages for 3D Mapping Example
library(geojsonio)
library(mapdeck)

Set your working directory

For this example, we will re-import all of the data and keep all of the species points to make our interactive map.The next few cells will be similar to what you already completed in the earlier exercise.

# Set the working/base directory where our project files are stored. This will be unique to your computer.
baseDir <- getwd()

# Create an output directory for any processed files
dir.create(path = "output")
outDir <- paste0(baseDir,"/output")

Import and set up the data

Again, we’ve done this, so we will do it in a single cell.

  1. Import the species sample location point data. /n
  2. Import the protected areas raster (raster) /n
  3. Convert the point data to a spatial points data frame (sp).
  4. Crop the raster to the extent of our point data (raster)
# 1. Read in the species sample points again...
d1 <- read.csv(paste0(baseDir,"/data/cucurbitaData.csv"))

# 2. Read in the protected areas raster again...
protLands <- raster::raster(x = paste0(baseDir,"/data/proAreas.tif"))

# 3. Convert the .csv to a spatial points data frame, again...
coordinates <- d1[,c(4,3)]
proj4 <- protLands@crs
sp1 <- sp::SpatialPointsDataFrame(coords = coordinates,
                                  data = d1,
                                  proj4string = proj4)

# 4. Crop raster to spatial data frame extent
protCrop <- raster::crop(x = protLands, y = sp1)

Extract protected area status to new column

Next we are going to extract information from the protected areas raster to a new column in our data frame called “protected”. It will fill the each row with a “1” if it falls within a protected area, and an “NA” if it falls outside of one. 

sp1$protected <- raster::extract(protCrop, sp1)
head(sp1)
##   X             taxon latitude longitude type protected
## 1 1 Cucurbita_cordata 28.95470 -113.5625    G         1
## 2 2 Cucurbita_cordata 24.28577 -111.0110    H        NA
## 3 3 Cucurbita_cordata 25.93470 -111.5421    H         1
## 4 4 Cucurbita_cordata 25.94749 -111.7300    H         1
## 5 5 Cucurbita_cordata 25.95765 -111.4609    H         1
## 6 6 Cucurbita_cordata 26.20671 -111.6175    H         1

Label all points based on whether they do, or don’t, fall in a protected area

Now we want to label points that fall within protected areas so we can map protected species locations. Since we already extracted protected area information to our data frame, records that don’t fall within protected areas will have an “NA” value. In this example we are going to change the “NA” values to “0”. Another option would be to simply omit all “NA” values from the table if we only wanted the protected points, but in this case we want to visualize the differences. 

sp1.class <- sp1

indexNA <- which(is.na(sp1.class$protected))

sp1.class[indexNA, "protected"] <- 0 
head(sp1.class)
##   X             taxon latitude longitude type protected
## 1 1 Cucurbita_cordata 28.95470 -113.5625    G         1
## 2 2 Cucurbita_cordata 24.28577 -111.0110    H         0
## 3 3 Cucurbita_cordata 25.93470 -111.5421    H         1
## 4 4 Cucurbita_cordata 25.94749 -111.7300    H         1
## 5 5 Cucurbita_cordata 25.95765 -111.4609    H         1
## 6 6 Cucurbita_cordata 26.20671 -111.6175    H         1
# We could also omit ALL "NA" values for an "all protected" data frame. But we can ignore this code for now.

# sp1.nona@data <- na.omit(sp1.nona@data)
# head(sp1.nona)

Load species rasters

Next we are going to load the species range maps, or species distribution model rasters found in your “/data” folder. We didn’t use these files in the previous example, but they are useful for mapping species ranges along with the sample point data. The file is in “.rda” format, so we will use the load() function and then pull each raster out of the file as individual raster objects. 

# Load the .rda file
path <- paste0(baseDir,"/unzip/CucurbitaRasters.rda")
load(path)

# Extract each species' raster to individual objects
cordata <- CucurbitaRasters$cordata
palmata <- CucurbitaRasters$palmata
digitata <- CucurbitaRasters$digitata

# Plot an example...
plot(digitata)

Filter raster data

We can see from the map above that the data contains presence and absence values for the range map. We only want to plot the presence data (i.e., “1”) in our web map, so we need to filter the data. 

digi.filter <- digitata
digi.filter[digi.filter[] < 1] = NA
digi.filter
## class      : RasterLayer 
## dimensions : 144, 150, 21600  (nrow, ncol, ncell)
## resolution : 0.1666667, 0.1666667  (x, y)
## extent     : -124.5, -99.5, 20.33333, 44.33333  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## source     : memory
## names      : digitata 
## values     : 1, 1  (min, max)
# Check the output...
plot(digi.filter)


Let’s do that for the other species range maps…

corda.filter <- cordata
corda.filter[corda.filter[] < 1] = NA

palma.filter <- palmata
palma.filter[palma.filter[] < 1] = NA

Create a Leaflet Web Application

The previous tutorial had us create a simple Leaflet web map as an example, but Leaflet has so much to offer! The following cell will create a fully functional web application, where we can view both point vector and raster layers, make a heatmap of location density, and visualize different base maps. This section is meant to provide a deeper overview of the capabilities of Leaflet for making interactive web maps and applications. Again, one exciting, but very easy feature to add is the ability to easily make heatmaps based on the point location densities.

Choosing colors with Color Brewer

A great place to choose color values (HEX, RGB, CMYK), especially for finding colorblind safe and print optimized color values, is from is Color Brewer

Leaflet Provider Basemap Tiles

Here we are grabbing Leaflet provider basemap tiles. You can see all of the options here.

Other map options…

There are many (though not unlimited) options for customizing your web map by utilizing all the features in the Leaflet package. You can extend the map’s functionality with a growing ecosystem of plugins and add-on packages (e.g., leaflet.extras), or just join a future tutorial on JavaScript. 

pal <- colorFactor(c("#5ab4ac", "#d95f02", "#7570b3"), domain = c("Cucurbita_cordata", "Cucurbita_palmata", "Cucurbita_digitata"))
propal <- colorFactor(c("#e31a1c", "#1f78b4"), domain = c(0, 1))

map <- 
  sp1.class %>% 
  
  leaflet() %>%
  
  addProviderTiles(providers$CartoDB.Positron, group = "CartoDB", options = providerTileOptions(minZoom = 4, maxZoom = 12)) %>%
  addProviderTiles(providers$OpenTopoMap, group = "Topographic", options = providerTileOptions(minZoom = 4, maxZoom = 12, opacity = 0.45)) %>%
  addProviderTiles(providers$Esri.WorldImagery, group = "Satellite", options = providerTileOptions(minZoom = 4, maxZoom = 12)) %>%

  addCircleMarkers(color = ~pal(taxon), radius = 4, weight = 1, fillOpacity = 0.60, group = "Species Records") %>%
  addCircleMarkers(color = ~propal(protected), radius = 4, weight = 1, fillOpacity = 0.60, group = "Protected Status") %>%

  addHeatmap(lng=~as.numeric(longitude),
             lat=~as.numeric(latitude),
             radius = 8,
             group = "Density (Heatmap)") %>%
  
  addRasterImage(protCrop, colors = "#fdbf6f", opacity = 0.75, group = "Protected Areas") %>%
  addRasterImage(corda.filter, colors = "#1b9e77", opacity = 0.65, group = "Cordata Range") %>%
  addRasterImage(digi.filter, colors = "#d95f02", opacity = 0.65, group = "Digitata Range") %>%
  addRasterImage(palma.filter, colors = "#7570b3", opacity = 0.65, group = "Palmata Range") %>%
  
  addLegend("bottomleft", pal = pal, values = ~taxon,
    title = "Taxonomic Name",
    opacity = 1, group = "Species Records") %>%
  
  addLegend("bottomleft", pal = propal, values = ~protected,
    title = "Protected Status",
    opacity = 1, group = "Protected Status") %>%

  addLayersControl(
    baseGroups = c("CartoDB (default)", "Topographic", "Satellite"),
    overlayGroups = c("Species Records", "Density (Heatmap)", "Protected Status", "Protected Areas", "Cordata Range", "Digitata Range", "Palmata Range"),
    options = layersControlOptions(collapsed = TRUE)) %>%
  
  # Set the initial view and zoom level of the map
  setView(-115, 30.75, zoom = 5) %>%
  
  # Using hideGroup gives us the option to hid map layers on launch to declutter our map
  hideGroup("Density (Heatmap)") %>% 
  hideGroup("Protected Status") %>% 
  hideGroup("Protected Areas") %>% 
  hideGroup("Cordata Range") %>% 
  hideGroup("Digitata Range") %>% 
  hideGroup("Palmata Range") %>%
  
  # This button add on allows us to reset the map to the original extent
  addResetMapButton()
  
map

Save the web map as an .html file

Finally, we can save our map as an .html file so we can explore the map in a web browser!

saveWidget(map, file= paste0(outDir,"leafletMap.html"))

Check your “/outputs” folder for the .html file. Double click and see how the map looks in your web browser of choice.

Create a 3D Web Map

Okay, one more fun map experiment! Let’s use the mapdeck package to make a 3D web visualization of our data. This package uses “Deck.gl”, a user-friendly WebGL javascript library that integrates nicely with Mapbox, which, like Leaflet, is a great web mapping platform.

WebGL is a JavaScript API for rendering interactive 2D and 3D graphics in your web browser. Make sure your browser is WebGL compatible

WebGL Platforms worth exploring:

Extract spatial data to polygon

Here I want to extract the species sample points to ecoregions so that we have something informative to plot. It’s completely possible to do this with our available data in R, but we only have 15 minutes, so I’ve already prepared a shapefile in ArcGIS Pro with species location attributes added. We haven’t loaded a shapefile into R yet, so let’s do that and quickly view the data. 

# Import and view 
ecoCount <- rgdal::readOGR(paste0(outDir,"/eco_feature_count_04162021.shp"))
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\niko8816\Desktop\NK\Jobs\CSU_Geospatial_Centroid\Research-and-Program-Coordinator-main\output\eco_feature_count_04162021.shp", layer: "eco_feature_count_04162021"
## with 18 features
## It has 4 fields
plot(ecoCount, 
     col=ecoCount$COUNTS,
     main = "Species Sample Count By Ecoregion")

Now we need to convert the shapefile to a format that mapdeck can read.

# Convert shapefile to GeoJSON
ecoCount.json <- geojson_json(ecoCount)

ecoCount.sf <- geojson_sf(ecoCount.json)
ecoCount.sf
## Simple feature collection with 18 features and 4 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -123.1593 ymin: 20.41155 xmax: -99.53937 ymax: 42.06112
## Geodetic CRS:  WGS 84
## First 10 features:
##                          ECO_NAME COUNTS latitude longitude
## 1                Apache Highlands     88 31.56698 -111.0514
## 2    Arizona-New Mexico Mountains      2 32.96590 -108.5767
## 3          Baja California Desert     96 24.28577 -111.0110
## 4        California Central Coast     13 35.60691 -120.1635
## 5          California South Coast    340 31.66667 -115.9500
## 6               Chihuahuan Desert     31 32.23387 -108.0502
## 7                Colorado Plateau      1 35.82304 -113.5702
## 8                     Great Basin      1 38.10881 -119.0518
## 9            Great Central Valley     62 35.21665 -118.7277
## 10 Gulf Of California Xeric Scrub     38 28.95470 -113.5625
##                          geometry
## 1  MULTIPOLYGON (((-109.0889 2...
## 2  MULTIPOLYGON (((-105.8239 3...
## 3  MULTIPOLYGON (((-111.6571 2...
## 4  MULTIPOLYGON (((-120.2269 3...
## 5  MULTIPOLYGON (((-119.373 34...
## 6  MULTIPOLYGON (((-99.57135 2...
## 7  MULTIPOLYGON (((-112.3335 3...
## 8  MULTIPOLYGON (((-113.1648 3...
## 9  MULTIPOLYGON (((-118.7917 3...
## 10 MULTIPOLYGON (((-110.3555 2...

3D Web Map

Now let’s plot some data to see how it looks. First we will test just using a simple data frame of species sample locations.

This package uses the Mapbox API for the basemap, so we will need to register for a free Mapbox account to access our user token. I’ve included my public token for this example, but please be kind and create your own for any future exploration. 

# Set the Mapbox key
key <- set_token("pk.eyJ1Ijoia290bGluaWMiLCJhIjoiZ3FCVFdTZyJ9.WfWEgjDKl9N0NKUyeLlAgA")

mapdeck(
    style = mapdeck_style('dark'),
    pitch = 45,
    location =  c(-112, 31),
    zoom = 5) %>%
  
  add_grid(
    data = d1,
    lat = "latitude",
    lon = "longitude",
    cell_size = 5000,
    elevation_scale = 200,
    layer_id = "taxon",
    update_view = FALSE
  ) %>%
  
  add_heatmap(
    data = d1,
    lat = "latitude",
    lon = "longitude",
    weight = "X",
    colour_range = colourvalues::colour_values(1:6, palette = "viridis"),
    update_view = FALSE
  )
# We are going to exaggerate our data values so that the features "pop" more. 
ecoCount.sf$COUNTS <- ecoCount.sf$COUNTS * 100
mapdeck(
    style = mapdeck_style('dark'),
    pitch = 45,
    location =  c(-112, 31)
  ) %>%
  
    add_polygon(
    data = ecoCount.sf,
    layer = "polygon_layer",
    fill_colour = "ECO_NAME",
    elevation = "COUNTS",
    legend = TRUE
  )

Now run the following cell to save the map as an .html document

ecodeck <- mapdeck(token = key, style = mapdeck_style("dark")) %>%
    add_polygon(
    data = ecoCount.sf,
    layer = "polygon_layer",
    fill_colour = "ECO_NAME",
    elevation = "COUNTS",
    legend = TRUE
  )

# We are going to save the widget as an .html file again so we can open it in our browser.
saveWidget(ecodeck, file="3Dmap.html")

Congratulations

Good work! You’ve finished this additional portion of the tutorial.